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Abstract 

We consider the nonlinear extension of the Kuramoto model of globally coupled 
phase oscillators where the phase shift in the coupling function depends on the 
order parameter. A bifurcation analysis of the transition from fully synchronous 
state to partial synchrony is performed. We demonstrate that for small ensembles 
it is typically mediated by stable cluster states, that disappear with creation of 
heteroclinic cycles, while for a larger number of oscillators a direct transition from 
full synchrony to a periodic or quasiperiodic regime occurs. 
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1 Introduction 



A model of coupled limit cycle oscillators explains a variety of natural phenom- 
ena in various fields of science. The applications range from the description 
of the collective dynamics of Josephson junctions [TJ, lasers [2], and electro- 
chemical oscillators [3 J to neuronal populations jl], etc. Very often, when the 
oscillator network is not too sparse, it can be approximately considered as 
fully connected, or globally coupled. 
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Ensembles of weakly interacting units are successfully treated within the 
framework of phase approximation |5,6,7J. Most popular is the Kuramoto 
model of sine-coupled phase oscillators, or its extension, the Kuramoto-Saka- 
guchi model [SJ. This model explains self-synchronization and appearance of 
a collective mode (mean field) in an ensemble of generally non-identical el- 
ements; the transition to synchrony occurs at a certain critical value of the 
coupling constant that is roughly proportional to the width of the distribution 
of natural frequencies Pf6|9|T0] . 

An extension of the Kuramoto model for the case of nonlinear coupling has 
been suggested in our recent publications [TTfTS] . see also [T3]rhf] . Nonlinearity 
in this context means that the effect of the collective mode on an individual 
unit depends on the amplitude of this mode, so that, e.g., the interaction of the 
field and of a unit can be attractive for a weak field and repulsive for a strong 
one. Formally, this is represented by the dependence of the parameters of the 
Kuramoto-Sakaguchi model (the coupling strength and the phase shift) on the 
mean field amplitude. The model exhibits nontrivial effects like a destruction 
of a completely synchronous state and appearance of partial synchrony in 
an ensemble of identical units. Moreover, in this setup the frequencies of the 
collective mode and of oscillators can be different and incommensurate. 

An analytical description of the dynamics of oscillator ensembles remains an 
important and challenging problem. A seminal work in this direction is that 
of Watanabe and Strogatz (WS) [ToTTB] . The WS theory is a powerful tool 
that provides a nearly full dynamical description of ensembles of identical os- 
cillators, sine-coupled to a common external force. In particular, this force 
can be the mean field of the population, so that for the case of identical units 
the WS theory almost completely describes the Kuramoto-Sakaguchi and the 
nonlinear models (see [T2] ). This description is given in terms of three collec- 
tive (macroscopic) variables, hereafter called the WS variables, plus constants 
of motion. The collective variables obey 3 WS equations (see [IT]); thus, the 
dynamics of an ensemble of identical elements is effectively 3-dimensional. 
However, the WS theory has one drawback: it cannot describe certain cluster 
states, i.e. regimes where the oscillators build identical groups. In this paper 
we complement the WS theory by performing a direct bifurcation analysis of 
the dynamical phase equations of the model of nonlinearly coupled phase oscil- 
lators. We will especially emphasis on cluster states and their bifurcations, in 
particular on the heteroclinic cycles (see [18] for a recent review of robust het- 
eroclinic cycles) that can be hardly treated within the WS approach. We will 
see that the role of clusters is mostly important for small ensembles. Because 
of identity of the oscillators, the system possesses a permutation symmetry, 
so we employ the corresponding bifurcation approach (see, e.g., [T9f20] ). 

The paper is organized as follows. We introduce the basic model in Section |2j 
Then in Section[3]we discuss general properties of bifurcations, possible attrac- 
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tors and their interpretation as different synchronization patterns. In Section H] 
we present bifurcation diagrams for a model of nonlinearly coupled oscillators 
with quadratic nonlinearity [12] . In Conclusion a relation to the WS theory is 
discussed. 



2 Model of nonlinearly coupled phase oscillators 

We consider an ensemble of N limit cycle oscillators, described by their phases 
6i e [0, 2ir), i = 1,...,N. They are assumed to interact globally, via the 
complex mean field 

1 N 

having amplitude r and phase ip: 

9 i = u i + G(r,i(>,9 i ). (2) 

Here Ui are natural frequencies of the oscillators and G is the coupling func- 
tion. Different popular models correspond to different choices of coupling func- 
tion G. The case G(r,ip,9i) = rKlm(e l ^~ e ^) corresponds to the famous Ku- 
ramoto model [6], while the choice G(r,-ip,9i) = rKlm(e~ ia e l ^~ e ^) yields the 
Kuramoto-Sakaguchi model [8]. 

In this paper we focus on a coupling function that nonlinearly depends on the 
amplitude of the mean field r and on a set of parameters (3: 

G(r, ip, 9i) = rK(r, f3) sin(^ -9 + a(r, (3)) . (3) 

This model has been introduced in |TT] and studied in the thermodynamic 
limit iV — y oo in [12]. In this paper we focus on the properties of small en- 
sembles of nonlinearly coupled oscillators, restricting our analysis to the case 
of identical oscillators Ui = u and of phase nonlinearity only K(r, 0) = 1. 
The latter restriction is not very important, as the cases where K(r, (3) can 
change sign are in fact trivial. Substituting ([3]) in ([2]) we obtain an equivalent 
formulation of the ensemble dynamics 

1 N 

Oi = uj + —Y J MO j + a(r, (3)) . (4) 

ly 3=1 

To exploit the phase-shift symmetry of this system one can describe the system 
dynamics in terms of the phase differences 

<p i = 9 1 -9 i+1 , i = l,...,N-l, (5) 
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thus reducing this A^-dimensional system to the (N — l)-dimensional system 



1 



N-l 



wn(<pi - (pj + a(r, 

j=i,3& 

N-l 

sin(v9i + a(r,p))+ sin (^i ~ «( r > P)) 

7 = 1 



(6) 



One can check that order parameter r can be written in phase differences as 



1 

N 



\ 



N-l 



N + 2 Y (cos((fj) + cos((fi - ipj)) 



(7) 



Below we will discuss synchronization transitions in the system studying in- 
variant manifolds, fixed points, cycles, heteroclinic cycles and their bifurca- 
tions for the system in phase differences (jBj). 



Before proceeding to the analysis, we mention that system (J3J possesses sym- 
metries given by all permutations of the oscillators [21]. Due to identity of the 
oscillators, the main dynamical regimes appear as invariant sets of the system: 
1) Completely synchronous solution, where all the oscillators are in the same 
state: 



= {{9 U ---M : 0i = 2 
3) Completely asynchronous solution 



On}', 



N 



M = <(e x ,...,e N ): = 

3=1 



The set M. is a union of invariant manifolds of dimension N — 2 for N > 3 
[22] . it corresponds to the case of vanishing order parameter r = 0. 
2) Cluster states, where groups of oscillators have identical phases. A general 
n-cluster state can be written as (up to permutation of indices) 



N ) : 6>i 



Opi i Qpi+i 



7 Pl+P2 ) 



: Q-r—vn— 1 . -, 



(9) 

where p\ + p 2 + • ■ ■ + p n = N. We will be mainly interested in 2-cluster states 
(we will see that only such states appear as stationary solutions) 



V 2 = {(9 1 ,...,9 



N) ■ VI 



N } 



(10) 



characterized by the partition (p : N — p). 



3 General analysis of synchronization and bifurcations 

In this section we study general bifurcation scenarios in the system of nonlin- 
early coupled oscillators @, to be illustrated by particular examples in the 
next section. 

3.1 Bifurcations in the Kuramoto-Sakaguchi model 

We start with the simplest case of linearly coupled oscillators. Here model PJ 
reduces to the standard Kuramoto-Sakaguchi model which we write as 

N 

N9 t = gi (6 Xl ...,6 N ,a) = -J2 sin(0i - Q i ~ «) ( n ) 

i=i 

Equilibria. To describe the steady states of the corresponding system in 
differences ipi = Q\ — 1+i , we need to solves the system of iV — 1 algebraic 
equations 

gi(9 1 ,...,e N ,a)-g i (e 1 ,...,9 N ,a) = 0, i = 2,---,N, (12) 

where a is a scalar parameter. The next lemma helps us to characterize the 
steady states of the system (TTTT) . 

Lemma 1 For any a G T 1 , the set (6\, • • • , 9n) satisfies system of equations 
KW\) if and only if one of the following three conditions is fulfilled: 

1) 0\ = - — = On, 

2) Ef=i e** = 0, 

s) e l = e 2 = ■ ■ ■ = e p ^ e p+1 = e p+2 = ■ ■ ■ = e N , p = 1, • • • , n - 1, (plus ail 

possible permutations). 

This Lemma means that the only possible steady states are that of complete 
synchrony (one cluster), complete asynchrony, and of two clusters. 

Proof: It is easy to check that states 1) - 3) satisfy the system (\V2\i . We will 
show that the roots of the system (1121) satisfy 1) - 3). We can re- write (1121) in 
the following way: 

N N 

(sin(#i — a) — sin(#j — a)) cos 9j — (cos(6 l i — a) — cos(6'i — a)) sin 8j = 0, 

i=i j=i 

(13) 

where i = 2, . . . , N. We consider four possible cases. 

A. If Y,jLi si 11 9j — and Y,f=i cos Qj — simultaneously, then the condition 
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2) of the lemma satisfies. 

B. The next possible case is that of Y^jLi^Oj — but J2f=i cos ^j 7^ 0. In 
this case (TTB"j) implies 

sin(#i — a) = sin(0j — a), i = 2, . . . , N. 

The last system shows that we can obtain only two-cluster solutions: 




-0i + 2a + vr, i=p+l,...,N, 



which must satisfy equations 

p sin 9i - (N - p) sin(0i - 2a) = 0, p = 1, . . . , N. 

The last equations arise from Y^f=i sin 9j = and they show that two-cluster 
states are possible only for some values of parameter a in this case. Note that 
the case p = N corresponds to a one-cluster solution (condition 3 reduces to 
condition 1). 

C. Consider the case, where X/^=i c °s0j — 0, and X}jLi sm 0j ^ 0. As in the 
previous case we obtain the possibility of two-cluster (or one-cluster, if p = N) 
states only: 

{9 1 , i = 2,...,p, 

-01 + 2a, i=p+l,...,N, 

which satisfy conditions 

p cos 0i + (N - p) cos(0i - 2a) = 0, p = 1, . . . , N. 

D. Consider (0i, • • • , 9 n ) such that Y,f=i sinOj ^ and J2j=x cos 0j ^ 0. Denote 
S := X/Jli sin0j, C := X/Jli cos0,,-, s" := sin(0j — a), cj := cos(0j — a). Then 
equation (FIB1) has the following form: 

K-OC-(c?-O5 = 0, i = 2,... (14) 

Dl. Suppose that sf — sf = for alH = 2, • ■ ■ ,N. Then using inequality 5^0 
we obtain c" — cf = for i = 2, • • ■ ,N. Equalities for and cj considered 
together yield e*^ j " a ^ — e*^ 1- ") = 0, j = 2, . . . , N, what means that all the 
values of 6j, j — 1, . . . , N, are equal. 

D2. Now let us consider another case, when there exists a number i such that 
s i — s io 7^ 0- Without loss of generality we can set zo = 2. Then from the first 
of equations (JHJ) we obtain 

c = K-c«)s/K-4). 
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Substituting C into the second equation of (|T4j) . we get 

S(s« - s «)(<£ - c«)/{s« - 4) - S{<% - c«) = 0. 
Using conditions S ^ and (s" — s 2 ) 7^ 0, we obtain 

K - «3)(c? - <£) - (sf - 4)(c? - 4) = 0, 

and then 

(«i - «) + - c£a?) + ( a «c£ - c«s%) = 0. 

After returning to the old notations and some transformations, we obtain the 
expression 

sin(0i - 9 3 ) + sm(6 2 - 9 l ) + sin(# 3 - 9 2 ) = 0, 

which already does not contain parameter a. We provide the last part of the 
proof by contradiction. The case D supposes that condition 2) is not valid. 
Now suppose that the conditions 1) and 3) are not satisfied as well. This means 
that there exists a solution (#1, . . . , 9n) of the system (fl3l) such that at least 
three variables 9 ix , 9 i2 , 9 i3 of this solution are not equal to each other. Without 
loss of generality we can set ii — 1, i 2 — 2, i 3 — 3 because we can replace 
variables using permutation (network has Sn symmetry). Inequalities 9\ 7^ 9 2 , 
9i ^ 3 , 9 2 ^ 9 3 imply that 

Bin(0! - 9 3 ) + sm(9 3 - 9 2 ) + sin(0 2 - 9 X ) = 

. . ( e 1 -e 3 \ . ( e 3 -e 2 \ . ( e 2 -e 1 \ 

= -4sm j sm j sm j ? 

This contradiction proves validity of either 1) or 3). 

D3. Consider a situation, when c" — c" o 7^ for some number i . In the same 
way as in the previous case D2 we prove that solutions of f|T3|) satisfy one of 
the conditions 1) or 3). 
Lemma is proved. 



Corollaries of Lemma 1. Lemma 1 implies that all steady states of the 
Kuramoto-Sakaguchi system, in terms of the phase differences, are one-cluster, 
two-cluster, or completely desynchronized states. As the two-cluster states 
constitute straight lines (plus those obtained by permutations of the variables) 

fl = <f2 = ■ ■ ■ = <fp 7^ Pp+l = Wf2 = • • • = lf N -! = 0, p = 1, . . . , N - 1, 

(15) 

all bifurcations of cluster steady states in this case are one-dimensional (in 
the sense that the normal forms are one-dimensional). Furthermore, to study 
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the existence of nontrivial cluster steady states we only need to solve scalar 
algebraic equations 



p sm(<p k — a) + (N — p) sin(<y9 fc + a) — (N — 2p) sin a = 0, k = 1, . . . ,p. (16) 
This equation has only two solutions on T 1 : ip k = and 



arppn „ | 2p(N-p)+(N 2 -2p(N-p)) cos(2 a ) ' 
cllCCUb 1 N 2 +2p(N-p)(cos(2a)-l) ' 



0;f ] U 



7T — 

' 2 



= (17) 



arccos " 1 ' 



/ 2p(N-p)+(N 2 -2p(N-p))cos(2a) \ U. \ , . 

V Af 2 +2p(AT-p)(cos(2a)-l) J " fc [ 2 ' V 



f;27r 



We can see that a bifurcation in the system ( 1111) occurs only when a = tt/2 
and it is transcritical. The bifurcation value of parameter a doesn't depend 
on the number of oscillators iV or on the cluster partition (number p). 

Note that in the case of a symmetric partition iV = 2p, equations fflBT) have a 
very simple form 

2psimp k cos a = 0. 

For this partition the only steady states are ip k = or ip k = vr, provided 
a 7^ ±7r/2. There is no any bifurcation on these lines for these values of the 
parameter. Vice verse, for a = ±7r/2 the whole two-cluster invariant line in 
the case of symmetric partition consists of fixed points. These fixed points are 
degenerate saddles (in the direction of lines with symmetry mentioned) and 
together with their one-dimensional manifolds they build a set of heteroclinic 
cycles. 

As it follows from lemma [1] and formula (fTT|) . the standard Kuramoto model 
of identical oscillators (a = 0) and the system with coupling a = ±ir/2 have 
a simple structure of the steady states. The standard Kuramoto model has 
only equilibria of two types: (i) equilibria that compose the manifold A4 (with 
vanishing order parameter) and (ii) equilibria that have coordinates differences 
Oj—9], equal to or to ±7r. In the case a = ±ir/2 all equilibria satisfy either (i) 
lie on the manifold Ai or (ii) correspond to a completely synchronous state, 
where 9j = 6^, j, k = 1, N, and r = 1. 



3.2 Bifurcations in model of nonlinearly coupled oscillators 



The model of our main interest (system (j3J) or, equivalently, (jSJ)) differs from 
the Kuramoto-Sakaguchi model only by the nontrivial phase shift a. Fortu- 
nately, using Lemma 1 we can localize steady states in a system of equations 
even more general than ( ITT]) , with r.h.s. containing an arbitrary scalar function 

a = a(0i, • • -,9 N ,f3), 
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where (3 is some vector of parameters (3 = . . . , (3 m ), m > 1. To do this we 
need to describe all solutions of the algebraic system 

g 1 (9 1 , . . . , 9 N , a(0i, ...,9 N , 0))-g i (9 1 , ...,6 N , a(0 1} . . . ,9 N , /?)) = 0, z = 2, ■ • • , N. 

(18) 

Lemma 2 ,0n) satisfy system (Tj2|) / or an V smooth scalar function 

a(9i, . . . , #at, /3) and vector of parameters (3 G M m z'/ and on/y if they satisfy 
one of the following conditions: 

1) 9\ — ■ ■ ■ — 9^, 

2) Ef=i e i9 i = 0, 

3) 0i = 2 = • • • = 9 P ^ 9 P+1 = 9 p+2 = ■ ■ ■ = 9 N , v = 1, ■ ■ • , N - 1, up to 
permutations. 

Proof: Let us assume that conditions of lemma [2] are violated for some 
fixed value of variables (9\, ■ ■ ■ ,9n) = • • • ,9%) and parameters f3 = (3° = 
(ft®, • • • , /3^J. Then lemma[T]is not valid for system fll2p for the fixed parameter 
value a = a (9i, • • • , 9%, (3°). This contradiction proves lemma [23 

Note that here we don't require from function a (and thus from coupling func- 
tion of the whole system) any type of symmetry. Nevertheless, all equilibria 
bifurcations are one-dimensional and they occur on the straight lines which 
are invariant for the system and are described by (fl6|) . However, in the paper 
we will consider coupling function g with permutation symmetry Sn and will 
describe bifurcations of the system using this symmetry property. 

Lemma [2] shows that (like in the standard Kuramoto-Sakaguchi model) all 
steady states of system (where a = a(r, (3)) belong only to the invariant 
manifold Ai or to clusters with isotropy S p x Sn- p . In the latter case the 
problem reduces to solving scalar algebraic equations 

psin(Lp k -a(r((p k ), f3))+(N-p) sin(v? fc +a(r(y9 fc ), (3))-(N-2p) sin(a(r(y2 fc ), (3)) = 

(19) 

for these steady states. In these equations the mean field amplitude r is defined 
according to (fT5l) and it depends only on one variable ip k , where k — 1, . . . , p. 
Also we can see that all steady state bifurcations have 1-dimensional normal 
forms for the model. Below we describe these and other bifurcations, illustrat- 
ing them with cases N = 3 and N = 4 (Figs. [T|2|) . 

Bifurcations of the completely synchronous state ipj = 0. The origin 
of the system is an equilibrium for any value of the function a(r,/3). Con- 
sider the Jacobian matrix of this system at the point (fj — 0, j — 1, • • • ,N — 1. 
All eigenvalues of this matrix have the same value: 

\ i = -Ncos(a(l,(3)), i = !,-•• ,N-1. 
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Fig. 1. Illustration of bifurcations of steady states for N = 3 (left panel in (a), the 
system in terms of phase differences (p is two-dimensional), and N = 4 (right panel 
in (a), the system in terms of phase differences cp is three-dimensional). Panel (b) 
illustrates particular transitions in the selected regions of the phase space (see text 
for details). 

This means that the origin of the system changes its stability when a(l, (3) = 
±7r/2. Also, as it was argued above, a bifurcation must be one-dimensional on 
each of the invariant lines with symmetry S p x SV-p- This bifurcation can be 
either a transcritical or a pitchfork one. A pitchfork bifurcation can happen 
only in the case of an even number of oscillators and this bifurcation occurs 
along invariant lines with the symmetry Sjv/2 x Sn/2 as it was shown in the 
work of Ashwin and Swift [2"T] . 

Thus, a typical bifurcation of the completely synchronous state is a transcrit- 
ical bifurcation (see raw (1) in Fig. [U(b)). These bifurcations occur simultane- 
ously on all invariant lines with the isotropies S p x Sn_ p , P 7^ N/2. Bifurcation 
parameters (3 = (/3i, . . . , (3 m ) are defined from the expression a(l, 0) = ±7r/2. 
The steady state at the bifurcation point is a degenerate saddle (all eigenvalues 
of the linearized system are zero). Y}j^ 1 ^ 2 1 C 3 N saddles, where [N] - integer 
part of N, meet together at the origin. The bifurcation changes stability of 
the origin along each of the one- dimensional directions. 

A pitchfork bifurcation of the origin (see raw (3) in Fig. []Jb)) occurs simul- 
taneously with the transcritical bifurcation, when the number of oscillators is 
even. Two saddles appear (disappear) from the origin (stable or unstable) and 
move in opposite directions along the lines which have Sn/2 x <Sjv/2 isotropy. 
In the case of an even N these saddles are usually generators of trajectories 
(one-dimensional manifolds) which can be parts of heteroclinic cycles, under 
some additional conditions. 



Clusters and their bifurcations. To find all other steady states on the 
invariant lines with the symmetries S p x SV-p we should solve appropriate 
algebraic system (TTS]) that satisfies (TH)]) - this means that we need to solve 
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Fig. 2. Illustration of bifurcations via heteroclinic cycles. Upper panel: HC appears 
via a saddle-node bifurcation and gives rise to a limit cycle. Middle panel: The case 
of Kuramoto-Sakaguchi, here at the bifurcation point a family of neutral cycles exist, 
while beyond it only a fully asynchronous steady state is stable. Bottom panel: HC 
appears via a transcritical bifurcation and gives rise to a limit cycle. 

one algebraic equation. Typically, steady states appear (or disappear) by pairs 
on each invariant lines for (fj e (0, 2n), j — 1, • • • , N — 1, and this appear- 
ance (disappearance) corresponds to a saddle-node bifurcation (see raw (2) in 
Fig. Hb)). A saddle-node bifurcation that occurs in the (N — l)-dimensional 
space (where our reduced system is considered) leads to the appearance of two 
new points, i.e. of two new two-cluster states. These two points have opposite 
stabilities along the one-dimensional manifold with isotropy S p x Sn- p , but 
the same stabilities transversal to these one-dimensional manifolds. In partic- 
ular, one of these two newly appeared points can be a stable or an unstable 
node. A stable node on the one-dimensional invariant line corresponds to a 
two-cluster with symmetry S p x Sn- p - 



Heteroclinic and limit cycles. Saddle steady states that appear in the 
transcritical and saddle-node bifurcations described above may have unstable 
manifolds that are connected to each other, thus constituting a heteroclinic 
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cycle (see also a similar structure described in [22]). When the heteroclinic 
cycle disappears, a usual limit cycle may appear corresponding to a periodic 
non-synchronous regime in system ($5$). We illustrate two types of such a bi- 
furcation in Fig. [2J In the upper panel we show an appearance of a limit cycle 
via a heteroclinic one, that appears at the saddle-node collision. Generation 
of a limit cycle by a saddle-node bifurcation via a heteroclinic cycle is a typ- 
ical situation in the system fl6]). In the case of an even number of oscillators, 
saddle-node bifurcations on the invariant lines can give possibility to connect 
different one-dimensional manifolds of saddles (pairs of saddles) generated by 
pitchfork bifurcation from the origin. The bottom panel illustrates a hetero- 
clinic cycle appearing at a transcritical bifurcation at the origin. Heteroclinic 
(or homoclinic) cycle consists of the origin point and loops of S p x Sn- p invari- 
ant lines. The middle panel in Fig. [2] shows the same transcritical bifurcation 
in the Kuramoto-Sakaguchi model. 

Another possibilities of a limit cycle to appear are an Andronov-Hopf bi- 
furcation of the point of M. (we will show this below), and a saddle-node 
bifurcation of two limit cycles. The existence of more complicated structures 
such as of a quasi-periodic torus or chaotic attractors is impossible for the 
system in phase differences ([6]). This follows from the Watanabe-Strogatz the- 
ory [16]. As it was shown in [T6lfT2] . the system (TjJ can be reduced to a skew 
three-dimensional system where the equation for one variable fully depends 
on two other ones. Thus, the dynamics of the two "driving" variables can be 
at most periodic, and the full dynamics at most quasiperiodic. In the terms of 
variables we use here, the "driving" variables correspond to phase differences 
tpk, their dynamics thus can be at most periodic. The full dynamics of phases 
9k includes one more integration and can be at most quasiperiodic. 



Multistability. If a saddle-node bifurcation generates a stable node while a 
stable node at the origin still exists, we obtain a bistability of a fully synchro- 
nized and two-cluster regimes. Note that depending on function a(r, (3), we can 
obtain many stable nodes on the invariant lines, resulting in a multistability 
of synchronous and different two-cluster states. 

One can also observe a coexistence of limit cycles appeared via different 
saddle-node bifurcations accompanied by heteroclinic cycles. Part of these 
cycles are stable but other ones are not. 



Attractors. As a result, one can observe the following types of possible 
stable regimes or their combinations in system ([6]): 

1) Complete synchrony ipj = 0, j = 1, • • • N. 
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2) Two-cluster regime with symmetry S p x Sn^ p . 

3) Limit cycle. 

4) Heteroclinic cycle. 

5) Manifold M (JV) . 



Stability of M.. Consider the invariant set Ai. This set is (N — 3)- 
dimensional in T^ -1 and consists of steady states of the system. To describe 
local bifurcations, we need to consider the property of Jacobin matrix 



d(gi(<pi, ■ ■ . ,y?jy-i,a), ■ ■ .,g N -i(<pi, . . .,<p N -i,a)) 
d((fi, . . .,<pn-i) 



on the points of the manifold At. We will show that N — 3 eigenvalues of 
Jacobian vanish, so there is no any motion inside the manifold. 

Lemma 3 Jacobian rank of the system is: 



rank( J) 



1, for 2-clusters with symmetry Sn/2 X SW/2, 

2, in other cases. 



Proof: Jacobian matrix J has the elements 



Jki 



Jkk — 



dgk _ _ 
dipt 



dgk 



dipk 
da 



da N ~ l 

cos(v9 fc - a) - - — V cos(v5 i - a) + 



\ I N-1 

- J cos(v9 fc + a) + ^(^fc - <fj + a) 

o W— 1 



da N ~ 1 

cos((pi - a) - cos(v3 fc -ifi + a) - — ^ cos(^- - a) + 



+ 



da 



I N-1 N 

cos(> fc + a) + Yl cos (fk ~ fj + «) 

\ j=l,j^k J _ 

Since we consider manifold Ai, then using (J8J) and (J5j) we obtain 

- a) = 



N-1 



cos a + V] cos((fj 

3=1 



and 



JV-l 

cos(<y9fc + a) + cos a + cos(y?fc — </?j + a) = 0. 

j=id¥=k 
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Thus in this case elements of Jacobian matrix are 



dg k 



dcpj 



= cos((fj -<pk- oc(0, (3)) - cos(v3j - a(0, /?)), j, k = 1, . . . , N - 1. 



Denote each column of matrix J by Jk, k — 1, . . . , N — 1. To prove that rank 
of matrix is not greater than two, we need to show that there exists the linear 
dependence between any three columns Jj, Jk, Ji of matrix J, i.e. there exist 
two scalar functions jj and 7^ such that 

IjJj + IkJk = Ji- 

One can check that the last expression is satisfied with functions 

_ sin(y?; - ip k ) _ sin(y?; - <pj) 

7j sin(v? j - ip k ) ' lh sin(v9 fc - ipj) ' 

when ipj 7^ ifk. Thus rank (J) < 2. 

We can rewrite equation for vectors in the form 

sm(ip k - (pi)Jj + sin((pi - ipj) J k + sin(^- - ip k )Ji = 0. 

All coefficients are not equal to zero in this expression when (fi ^ (p^ 7^ <pi. 
Thus rank (J) is not less than two when the system has at least a three-cluster 
regime. Then rank(J) = 2 for three-or-more cluster regimes. 

In the case of even number of oscillators N = 2p the system can have two- 
cluster states with symmetry Sn/ 2 x Sjv/2 that belong to invariant manifold 
M.. This means that in the last equation one coefficient is equal to zero that 
implies rank(J) = 1. The lemma is proved. 

The lemma shows that the Jacobian has iV — 3 eigenvalues equal to zero at 
the points of manifold M.. However, as it was shown, the rank of the Jacobian 
depends on values of the variables (i.e. on the coordinates of points on the 
manifold). Thus, to find the eigenvalues of J we need to consider not only 
2x2 minor of Jacobian matrix but the whole matrix. Each of the eigenvalues 
is a function of iV — 3 variables in the points of the manifold. Let us express the 
last two variables and <fN-i as the functions of the variables <pi, . . . , pns 
using expressions of real and imaginary parts of (jHJ). Then we obtain: 

(f N - 2 = arctan ( J -f ) - ^ arccos ( Jl - — - 1 j + ^(1 - sign(/i)), 




(fN-x = arctan ( ^ ) + ^ arccos ( J 1 - — - 1 j + — (1 - sign(/i)), 
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where 



N-3 N-3 

j=i i=i 



In the case of a uniform distribution of oscillators on the circle (splay state 
according to terminology used for Kuramoto model [16]), that is when 

^ = iV ' <Pj=j<Ph j = 2,...,N- 1, 



the eigenvalues of the Jacobian are 



Ajv-2,jv-i — ~2~( cos « ± « sm a). 

In general case the eigenvalues are 

N ( j 

Aiv-2,7v-i(v 9 i, • • -,<Pn-3) = y ( cosa ± ycos 2 a - h 2 (<pi, . . . , 

where • • • , </?at-3)| < 1 is some enough complicate smooth function. 

Therefore, we obtain an Andronov-Hopf bifurcation when a = ±7r/2. This 
bifurcation happens simultaneously in each point of manifold Ai except for 
the points with isotropy Sn/2 x Sjv/2 where function h(ip 1 , . . . , ^jv-3) = 0. 

In the case of three coupled phase oscillators, the zero-dimensional manifold 
consists of two points (2tt/3, 47t/3) and (47r/3, 2n/3). At the point of 
an Andronov-Hopf bifurcation each of these two points changes its stability 
and generates supercritically (or destroys subcritically) a limit cycle. With a 
further variation of a parameter, this limit cycle can grow in amplitude and 
disappear, either in a saddle-node/heteroclinic bifurcation, or via a saddle- 
node bifurcation of two limit cycles. 

More nontrivial situations can happen in the case of four globally coupled 
oscillators. Invariant manifold Ai^ in this case consists of six straight lines. 
Coordinates of such lines are (<pj, ir, Pj+ft) up to permutations. Invariant man- 
ifold has Z 2 isotropy. The function h(tpj) = sin ipj appears in the expressions 
for the eigenvalues of the Jacobian matrix. An Andronov-Hopf bifurcation 
happens simultaneously in each of the manifold points. Thus, we obtain a two 
dimensional surface that consists of limit cycles. Noteworthy, the bifurcation 
differs from a Neimark-Sacker bifurcation. Each of this cycles is attractive 
(repulsive) only inside the surfaces described by Watanabe-Strogatz theory 
[16], in other direction it is neutral. Possible way of this two-dimensional sur- 
face to disappear a is saddle-node heteroclinic bifurcation on invariant lines 
described. A two-dimension set of heteroclinic cycles occurs at the point of 
bifurcation (such a set was shown in figure 10 in |22j). 
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Another possibility is disappearance (appearance) of two limit cycles in a 
saddle-node bifurcation of cycles. We can note that such a bifurcation happens 
for each pair of limit cycles which belong to different two-dimensional sets of 
cycles. Such a bifurcation happens also inside the Watanabe-Strogatz surfaces. 
Thus we obtain a saddle-node bifurcation of two-dimensional surfaces, one of 
them is stable and other is unstable. 



4 Nonlinearly coupled oscillators with quadratic phase nonlinear- 
ity 

As an example of application of general picture outlined above we consider 
the model (jlj with particular dependence of the phase shift on the amplitude 
of the mean field [T2] : 

a = a{r,P) =/3 1 + /3 2 r 2 . (20) 

Here the two-dimensional space of parameters /3 2 ) is a cylinder IxTd 
(As> Pi), because the r.h.s. of the equations are 27r-periodic with respect to f3\. 
The oddness of the r.h.s. of the system implies the symmetry of the parameter 
plane (Pi,P 2 ) -»■ (-/^-AO- 
According to the consideration above, two types of bifurcation happen when 

cos(a) = cos(/?i + P 2 r 2 ) = 0. 

One possible bifurcation is an Andronov-Hopf bifurcation (AH) on the invari- 
ant manifold Ai. Since r = on the manifold, then we obtain two straight 
bifurcation lines 

Pi = tt/2 and pi = 3vr/2 (21) 

on the parameter cylinder. The line Pi = n/2 corresponds to a supercritical 
Andronov-Hopf bifurcation and the line pi = 3tt/2 corresponds to a sub- 
critical one. Manifold Ai is stable if pi G (7r/2, 3tt/2) and it is unstable if 
Pi E [0,7r/2) U (37r/2,27r). 

Another possible bifurcation is a transcritical bifurcation (TC) at the origin 
((Pj = 0, j = 1, . . . , N — 1) along each of the invariant lines with the symmetry 
S p x Sjy-p, P — 1, ■ ■ ■ , N — 1. A pitchfork bifurcation (PF) at the origin occurs, 
simultaneously with the transcritical one, along the invariant lines with the 
symmetry Sn/2 x Sn/2 for the system with even number of oscillators. At this 
bifurcation point the origin is a degenerate saddle with N — l zero eigenvalues. 
The origin point of the system corresponds to full synchronization state where 
order parameter r = 1. Therefore, straight lines 

p x + p 2 = tt/2 + vrm, me Z, (22) 
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correspond to the transcritical (TC) or to the transcritical-pitchfork (TC/PF) 
bifurcations in the parametric space. The origin (i.e. the regime of full syn- 
chrony) is stable when (3i + (3 2 G (— tt/2 + 27rm, ir/2 + 2irm), m G Z, and it 
is unstable when 0i + 02 £ ( 7r /2 + 27rm, 37r/2 + 2iTm). 

These two types of bifurcation are independent of the number of oscillators, 
so the grids of straight bifurcation lines (j21|22p are present at any bifurcation 
diagram model (120|) (Figs. [3], [T]) . Other bifurcations of the fixed points occur 
only on invariant lines ffl5|) and they all are of the saddle-node type. The 
expressions for the order parameter on the invariant lines (|15p are 



r 2 =r 2 (p,N-p) = —(2p(N-p) cos(^) + (N - l) 2 -2(p- l){N-p- 1) + 1), 



where k = 1, . . . is a variable that changes along invariant lines with 
x Sjv_ p isotropy. To find the coordinates of the steady states we need to solve 
equation ( fl9|) with this expression for the order parameter. This expression 
simplyfies in the case of even number oscillators to 



and is independent of the number of oscillators. Thus, it describes appearance 
(disappearance) of two points on the invariant line with symmetry Sn/ 2 X SV/2 
after each pitchfork bifurcation at the origin. The coordinates of these points 
(which are saddles) on the invariant lines are then 



These symmetric points are important because they are the basis for hetero- 
clinic cycles in cases of even number of oscillators. 

For any number of oscillators, together with the grid of the straight lines de- 
scribing Andronov-Hopf and transcritical bifurcations (I21|22p . there are a lot 
of bifurcation lines that correspond to saddle-node bifurcations (SN) on invari- 
ant lines. Some of these lines correspond to heteroclinic bifurcations (HC or 
HC(SN)), provided some additional conditions are satisfied. At these saddle- 
node/heteroclinic bifurcations limit cycles appear. Let us fix some parameter 
value Pi 7^ ±7r/2 and increase parameter 02- Then new and new saddle-node 
bifurcations occur on the invariant lines, steady states appear at these bifur- 
cations in such a way that their stability alternate along this line. Stability of 
appearing heteroclinic cycles also alternates with increasing of the parameter 
02- Thus, stabilities of limit cycles that move inside each invariant region after 
such a bifurcation also alternate. The period and the amplitude of each of the 
limit cycles decrease when parameter f3 2 increases, but no bifurcation of small 
limit cycle can happen because we demanded that 0\ was not equal to ±7r/2 
where such a bifurcation only could occur. 




k = 1,... 



■P 



(fk = i arccos — (tt(1 + 2m) — 2f3 1 — (3 2 ) 
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The number of the hyperbolic steady states increases with increasing of pa- 
rameter fa. These steady states tend to concentrate near the center part of 
the invariant lines. A saddle-node/heteroclinic bifurcation on the line with 
symmetry Si x SV-i usually occurs close to this central region, where the co- 
ordinate ip of the saddle-node point is close to tt. Thus, we can approximately 
calculate that such bifurcation occurs at 



The lines of a saddle-node/heteroclinic bifurcation alternates on the bifurca- 
tion cylinder. Each heteroclinic cycle generates a limit cycle (or a set of limit 
cycles for 4 and more oscillators) with the same stability. The period and the 
size of each cycle decrease with increasing of parameter fa. Stable and unsta- 
ble limit cycles enwrap each other inside invariant region, and alternate. A 
saddle-node bifurcation of limit cycles is impossible for the system considered 
because of monotonic increase of cycles sizes (it would be, however, possible 
for more complex than (120]) dependencies a(r), e.g. for a — fa + far 2 + far 4 ). 

Noteworthy, the first saddle-node bifurcation in the system usually happens 
when fa < 7r/2 — fa, and this bifurcation can generate a stable node. Then we 
obtain multistability of the fully synchronous state (the origin where r = 1) 
and the stable two-cluster states. If the first transcritical bifurcation (when 
fa > n / 2 — fa) doesn't produce a stable heteroclinic cycles, then the two- 
cluster states are the only attractors in the system. The stable nodes accumu- 
late on the invariant line with increasing fa, thus we obtain a multistability of 
two-cluster states when fa is large enough. The appearance of stable limit cy- 
cles after saddle-node/heteroclinic bifurcation eliminates one stable node (on 
each invariant line with the same symmetry). However, since a heteroclinic 
bifurcation happens more rarely than simple saddle-node bifurcations of two 
points, then the coexistence of a stable limit cycle with two-cluster states is 
typical for the system. Therefore, we can obtain multistability of all possible 
attractors in the system: full synchronous state, two-cluster state (with differ- 
ent order parameters), limit cycle of phase differences, heteroclinic cycle, and 
invariant manifold A4. 

4-1 Three interacting oscillators 

For three oscillators interacting according to phase shift (|2"U1) . the bifurcation 
diagram is depicted in Fig. [3J The corresponding bifurcations have been al- 
ready illustrated in Figs. [Tf2l above. According to this bifurcation diagram, we 
show in Fig. H] a schematic dependence of the synchronization states in the 
system as parameter fa changes while fa = 0. One can see that the basic tran- 
sition in terms of phase differences is: Full synchrony — > two-cluster state — > 
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Fig. 3. Bifurcation diagram for N = 3 oscillators in the parametric plane 

(cf. Fig. Ufa)). Right panel is an enlargement of the central part of the left one. TC 
- transcritical bifurcation (see Fig. dfb)(l)), SN — saddle-node bifurcation (see 
Fig. [Tfb)(2)), AH — supercritical Andronov-Hopf bifurcation (see Fig. [2(b)(4)), 
HC(SN), HC(TC) heteroclinic bifurcations (see Fig. [2}. Points A, B and C are 
co dimension-two bifurcation points. The region where a stable limit cycle exists 
(right panel) is surrounded by a supercritical AH bifurcation line and two lines of 
heteroclinic bifurcations of different types. 




Synchronous state Stable 

Manifold M Saddle 

Two-cluster state Unstable 
Limit cycle 



Fig. 4. Schematic bifurcation diagram for the (/32, r ) parametric plane and the case 
a = far 2 . TC — transcritical bifurcation, HC — heteroclinic cycle (bifurcation), 
SN — saddle-node bifurcation. 



periodic oscillations. The first transition is with hysteresis (i.e. in some small 
region of parameters full synchrony and two-cluster state coexist), and the 
second transition is via heteroclinic connection. On the diagram Fig. [3] several 
codimension-2 points are marked, we will discuss them in the next subsection. 
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4-2 Four and more coupled oscillators 



In the Figures Eland Owe show schematically a saddle-node/heteroclinic bifur- 
cation for the case N — 4. An unstable heteroclinic cycle (Figure [5]) consists 
of ten fixed points and ten 1-dimensional invariant manifolds that connect 
these points. Four saddle-node bifurcations happen simultaneously on four 
Si x S 3 invariant lines. This heteroclinic cycle includes also two saddles S' 
that belong to S2 x S2 invariant line. Thus, the system of four oscillators, that 
moves along this heteroclinic cycle, shows temporary switches between 1 + 3 
and 2 + 2 clustering. This unstable heteroclinic cycle is robust and it will exist 
also beyond the saddle-node bifurcation; however it will consist of two lines 
connecting S' only, like the stable cycle depicted in the figure. This stable 
heteroclinic cycle is shown inside the unstable one. It consists of two saddles 
S" and two connecting lines r 1; T 2 . The stable heteroclinic cycle appears at a 
saddle-node bifurcation on the invariant line in the same way as the unstable 
one, only for a smaller value of parameter fa. The next heteroclinic bifurcation 
will occur after merging of stable node N + and saddle S and it will produce 
causes stable heteroclinic cycle. 

Figure |6] shows an appearance of a 2-dimensional sets of stable limit cycles 
inside one invariant region. Four pairs of saddles S and stable nodes N (Fig. [6] 
a) ) collide and create 2-dimensional sets of heteroclinic cycles (Fig. |6]b) ). 
Beyond the bifurcation, when saddle-node points SN disappear, two hete- 
roclinic cycles appear, with 2-dimensional sets of limit cycles between them 
(Fig. [6]c) ). The set of limit cycles surrounds the 1-dimensional invariant set 
Ai. This set of limit cycles shrinks as parameter fa increases, but it never 
reaches manifold Ai. 

The heteroclinic cycles presented in Fig. lie on invariant surfaces and cor- 
respond to switches between the cluster states. They are borders of sets of 
limit cycles that exist inside the bulk of the phase space (that is bounded by 
the invariant lines and surfaces), and enwrap the manifold Ai. To an unstable 
HC corresponds a cylindrical set of unstable limit cycles, and to a stable HC 
corresponds a cylindrical set of stable limit cycles. In this way the structure 
of heteroclinic cycles determines the overall structure of the trajectories also 
outside of invariant manifolds. 

The bifurcation analysis of higher-dimensional cases (for N > 5) shows similar 
results. The bifurcation diagram consists of three types of lines: a straight line 
of an Andronov-Hopf bifurcation, a straight line of a transcritical bifurcation, 
and lines of saddle-node bifurcations on invariant lines with the symmetry 
S p x Sjv-p, p = 1, . . . , N — 1. All the saddle-node bifurcation lines have sim- 
ilar "tongue-like" form. The "tongue" is formed by two border lines: one lies 
left of the straight line of a transcritical bifurcation and approaches this line 
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Fig. 5. Schematic phase portrait for N = A coupled oscillators, showing only in- 
variant lines corresponding to cluster states 1 + 3 and 2 + 2, and planes connecting 
them. We illustrate a coexistence of two types of heteroclinic cycles. One unstable 
HC (blue) is shown at the bifurcation point, another stable HC (red) is beyond its 
bifurcation. 




a) b) 




Fig. 6. Schematic diagrams of the saddle-node/heteroclinic bifurcation for 4 globally 
coupled oscillators, a) - prior bifurcation, b) - at the bifurcation point, c) beyond 
bifurcation. 
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asymptotically with increasing of parameter fa, the second border line crosses 
the TC-lines. The saddle-node bifurcation generates a sink and source only 
on the invariant line with symmetry Si x <Sjv-i> while on other invariant lines 
with the symmetry S p x S^-p, V 1 a pah °f saddles appears. Therefore, 
there exist only 1 x (N — 1) stable clusters. Furthermore, heteroclinic cycles 
and stable limit cycles appear at a saddle-node bifurcation with Si x Sn-i 
symmetry only (the corresponding bifurcation lines are drawn with blue in 

Fig. ED- ' 

Let us discuss the codimension-two points marked in Fig. [3J At point C two 
borders of the saddle-node tongue meet. Only one stable state is involved in 
both saddle-node bifurcations (say, on the left line states 1 and 2 are created, 
while on the right line state 2 annihilates with state 3), at the codimension- 
two point C all three involved steady states meet. At another codimension-two 
point A the type of the saddle-node bifurcation changes. On one side (left to 
the point A) no heteroclinic cycle appears at the saddle-node, while right 
to point A the transition can be saddle- node/heteroclinic, provided — n/2 < 
fa (A) < 7r/2 and the line has symmetry Si x Sn-i- We have checked that the 
latter condition holds for N = 3, . . . ,8 only and for N > 9 one has fa(A) < 
— 7r/2. Thus, for N >9 there is no saddle- node/heteroclinic transition. Finally, 
the point B on the bifurcation diagram Fig. corresponds to a degenerate 
situation depicted in the middle panel of Fig. [2j 

We show bifurcation diagrams on the planes [fa, fa) for four and five coupled 
oscillators in Fig. The structure of these diagrams is basically the same as for 
three oscillators Fig. EJ but with some quantitative changes. To clarify these 
changes we compare in Fig.[H]the basic saddle-node "tongues" for N — 3, ... ,7. 
One can see that with increase of N the tip shifts down and for a fixed fa « 
the saddle-node bifurcation can be observed for a small number of oscillators 
only. Thus, for a fixed fa ~ the loss of full synchrony with increase of fa 
occurs as direct transition from full synchrony to periodic oscillations via a 
transcritical bifurcation (bottom raw in Fig. [2]), and not via clustered states. 



5 Conclusion 

In this paper we have performed a detailed bifurcation analysis of the nonlin- 
ear generalization of the Sakaguchi-Kuramoto model of globally coupled phase 
oscillators. The main novelty in addition to the consideration in the frame- 
work of WS theory [12J is the characterization of cluster states that in terms 
of phase differences appear (via a transcritical, a pitchfork, or a saddle-node 
bifurcation) as steady states on invariant lines of the corresponding cluster 
configurations. At saddle-node bifurcations these steady states disappear via 
heteroclinic cycles. Remarkably, heteroclinic cycles in this model are not de- 
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Fig. 7. Bifurcation diagrams for (a) N = 4 and (b) N = 5 oscillators in the (/?2 ? /3i) 
parametric plane. In (b): SN\ — the line of saddle-node bifurcation on S4 x Si 
invariant line, SN2 — the same with the 53 x 52 invariant line. 




Fig. 8. The comparative diagram of saddle-node bifurcation lines on the (/?2,/3i) 
parametric plane for the cases 3, 4, 5, 6, and 7 oscillators. The tip of the "tongues" 
shifts down as N grows, and for N > 9 is below thee line of Andronov-Hopf bifur- 
cation f3i = —ir/2. 

stroyed but remain to exist (for other examples of heteroclinic cycles in ensem- 
bles of identical phase oscillators see [23,24 20] ) . This is related to the partial 
integr ability of the system resulting from the WS theory. According to WS, 
because the equations have A^ — 3 constants of motion, periodic orbits form the 
families of corresponding dimensions, the heteroclinic cycles form the limiting 
cases of these families, describing cycles that include nearly-clustered states. 

The analysis performed in this paper complemented the conclusion on the 
transition from full to partial synchrony in nonlinearly coupled oscillator 
ensembles, made in [12]. We have demonstrated that for small ensembles 
the transition is of the type "full synchrony" — > "cluster state" — > "peri- 
odic/quasiperiodic partially synchronous state" occurs, while for a large num- 
ber of oscillators a direct transition "full synchrony" — > "periodic/quasiperiodic 
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partially synchronous state" is typical. 
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